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'P' Abstract 

Oh! 

r-| ! A fast and simple Monte Carlo program is presented that simulates single 

Bremsstrahlung in Bhabha scattering, e^e~ —>■ e~^e~'j, without constraints 
on scattering angles. This allows the study of this process at arbitrarily 
small, or even vanishing, scattering angles. Experimental cuts can be im- 
posed on an event-by-event basis, allowing for detailed studies of the process 
as a limitation to beam lifetimes, or a luminosity-measuring device, in e~^e~ 
storage rings. As an application, we show that the easy introduction of a cut- 
off parameter, corresponding to the characteristic distance between particles 
in the e^ bunches, gives a reduced cross section that is in good agreement 
with observation. 
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Program Summary 

Title of program: BBBREM 

Program obtainable from: R. Kleiss, NIKHEF-H, RO. Box 41882, 1009 DB 
Amsterdam, The Netherlands, t30@nikhefh.nikhef.nl 

Computers: systems supporting standard FORTRAN 77 

Programming language used: FORTRAN 77 

Memory required: about 200kb 

number of bits per word: 32 

Subprograms used: none 

Number of lines in distributed program: 386 

Keywords: Bhabha scattering, Bremsstrahlung, radiative processes, forward 
scattering, collinear singularities, Monte Carlo simulation, experimental cuts, 
beam lifetimes, luminosity monitoring 

Nature of physical problem: Radiative Bhabha scattering, e~^e~ -^ e"'"e~7, 
has a very large cross section at small scattering angles, and plays various 
roles in existing and future e~^e~ colliders. It can be an important background 
to several two-photon scattering processes; it forms the major ingredient in 
the finite lifetime of colliding beams; and it is a possible process by which 
the luminosity can be measured, by observing electrons or photons emerging 
at zero scattering angle. Accurate knowledge of its cross section is therefore 
important. 

Method of solution: Due to the extremely singular structure of the matrix 
elements, and the possibility of complicated or unusual experimental cuts, a 
straightforward integration of the cross section over the allowed phase space is 
impractical. We therefore construct a Monte Carlo algorithm that generates 
(random) events in phase space, with a distribution that matches the actual 



cross section as closely as possible. These events are assigned a weight which 
corrects for discrepancies between the actual and the approximate matrix 
elements: the average value of the weight in a sample of generated events is 
the Monte Carlo estimate of the cross section. Since each generated event is a 
complete description of the momenta of the produced particles, any conceiv- 
able experimental cut can be implemented, in addition to the single a-priori 
constraint, namely, a minimum value for the energy of the Bremsstrahlung 
photon. This value can be set (to essentially arbitrarily small values) by the 
user. By setting the weight of events that fail a particular set of cuts to zero, 
one obtains the cross section under those cuts. 

Typical running time: about 185 /xsec per generated event on SUN SPARC 
10; about 40 fisec per generated event on IBM 9000. 

Unusual features of the program: none 



Long Write-Up 
1 Introduction 

The process of radiative Bhabha scattering, 

e~{pi) e+(gi) -> e'{p2) e+(g2) l{k) , (1) 

where we have indicated the various momenta, is expected to play an impor- 
tant role in several physical problems. In the first place, it is the main process 
by which electrons are lost upon the collisions of e^ beams, and dominates 
the beam lifetime of LEP |jl]]. Secondly, it may be an important background 
to processes where single photons or tt^'s are observed at small angles, such 
as at the DA$NE collider 0. One therefore intends to install tagging devices 
that can, in principle, observe electrons and photons emitted at zero scat- 
tering angle. Finally, this last device can also serve as a luminosity monitor 



due to the very large event rate. It is obvious that the various cross sections 
related to this process have to be known accurately. 

Various authors have already treated the process (|^), using a variety of 
approximations to the matrix element. Since scattering over very small an- 
gles completely dominates the process, Bremsstrahlung off the positron and 
off the electron are essentially independent, and we use only the two Feyn- 
man diagrams where the electron is seen to radiate the photon. On the other 
hand, when all the particles' trajectories are essentially coUinear, the elec- 
tron mass m cannot be neglected with impunity. The result of the CALKUL 
collaboration 0, where the leading correction terms of order th? / E"^ are in- 
cluded (where E is the beam energy), is also inapplicable since it assumes 
that the electron mass is negligible with respect to t, the momentum transfer 
of the positron. Of the early publications, the work of Altarelli et al. must be 
mentioned. They computed the total cross section [Q and various differen- 
tial distributions 0. Similar results were obtained by a number of Russian 
authors |^. With the possibility of actually measuring the process (|l]) at 
small or zero angle, these analytic, totally inclusive, results are no longer 
adequate, and the use of approximate matrix elements somewhat doubtful. 
In a version of the matrix element, appropriate to small-angle scattering, 
and correct up to truly negligible terms, was described in detail. In two 
recent publications, and Italian collaboration has presented results for cross 
sections obeying a variety of angular and energy cuts 0, and described the 
integration program used to obtain these ^ (the matrix element employed 
in has, in fact, been checked to be completely identical to that of 0|). 

In this paper we describe another approach to the study of (|ID, namely 
that of a full-fledged Monte Carlo simulation. The user generates any desired 
number of events, that is, sets of momenta of the three produced particles. 
Their distribution in phase space is matched as closely as possible to the 
exact distribution. The remaining factor is included in each event's weight. 
The average of the weights in the sample is the Monte Carlo estimate of the 
cross section. This is a very flexible approach since any experimental cut 
can be imposed by simply setting the weight of events that fall outside the 
cuts to zero. Moreover, since the momenta are explicitly given, the produced 
particles can be tracked by other simulation programs, either into a detector, 
or in the accelerator. A good check is provided by the fact that those cuts 
that are provided for in |^ must also yield the cross section calculated by 



that program, by completely different meansQ. 

2 The Monte Carlo 

We shall now describe our approach to the cross section of Eq.(|l]). We start 
by defining the necessary variables. We shall use the following invariants: 
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Sl = (P2 + ?2 
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p'^k^ , A2 = p'ikf. 



(2) 



Together with the azimuthal orientation of the event around the beam axis, 
these serve to completely specify all the momenta. In terms of these invari- 
ants, the correct form of the matrix element, squared and summed/averaged 
over the particle spins, is given by 
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(3) 



where e is the electron charge. The distribution of the generated events over 
the phase space shall be proportional the approximate matrix element 
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where Ci and C2 are defined later on. This simple form captures the essentials 
of the process dynamics, and also, owing to clever definitions of the Ci,2, leads 
to a simple result when integrated over the whole phase space. Denoting this 



^A small discrepancy was actually observed, due to the fact that in 
netic coupling constant a used was 1/(138.318) instead of 1/(137.036) 



the clectromag- 



integral (or approximate total cross section) by aapp, we then define the 
weight of a generated event by 

. . M 

weight = — dapp , (5) 

J"app 

for those events that satisfy the desired cuts, and zero otherwise. Note here, 
that only one cut is always present, namely a lower limit on the photon en- 
ergy. This lower limit is given by kE, where < k < 1, and serves to avoid 
the infrared singularity lurking at /c° = 0. 

We shall deal with two Lorentz frames: the laboratory frame, and the R- 
frame, that is, the rest frame of R^ = P2 + k^^. The most convenient picture 
of the process ([l|) is then the emission of a virtual photon by the positron 
(which is best described in the lab frame), followed by 'Compton' scattering 
of the virtual photon and the electron (which finds its easiest description in 
the i?-frame). We therefore have as 'natural' invariants the following set: 

• the components of the positron momentum transfer g^ = Qi — q2, 
namely t = q^q^, y = q'^/E, and (j), the azimuthal angle of q2 around 
the beam axis. These are all defined in the lab frame. 

• the scattering polar angle 6^ and azimuthal angle (p-y of the photon with 
respect to the incoming electron momentum. These are defined in the 
i?-frame, with c-y = cos^^. 

In terms of these variables, the phase space integration element can be written 

as 



0?$ = d'^P25{pl - ma) d'^q25{ql - m^) d'^k5{k'^) 5'^{pi + qi - P2 - q2 - k) 
2 2 

= — dy dt d(f) dc^ d(j)y , (6) 

where w"^ = R^R^ = sy + m?. The only nontrivial phase space bounds are 
those on t\ we have t,w-,Y, < t < tina-x-, with 



tniin = 2m' -2Eq'i-2 {m' -EqlY ~m\E-q'i 
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m^sy'^/tmm ■ (7) 



Now, y is typically very small, of order vn? jE'^] therefore, |tmax| can become 
extremely small. 

We may now write the invariants Ai 2 in terms of the chosen variables: 

A2 = syjl , 

Ai = s?/A(w^,m^t)(e + t;^)/(4t(7^) , 



A(w^,m^, t) 



O O \ O / 

w + m —t) — Aw m , 

(w^ + m^ - t)/\{w^, m^, t) - 1 
1 — c^ . 



Now, the approximate differential cross section is 



da. 



app 



The angular integrals are now trivial, and we find 

'2 + e\ C1C2 



dy dt dc. 



daa^V = 4" log 



e / \t\y \(w'^ , m'^ , t) 
we now choose Ci so as to simplify this expression: 
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this approximation is justified when |t| is very small. We then have 
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Now we do the t integral; owing to our choice of Ci, it gives the fairly simple 
result 
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sy^ \ m^ 

X{w'^,m'^,t) + sy — t' 



[^(tmax) - J{tmm)] dy 



^\{w'^, m^, t) — sy — t ) 
It is now time to choose C2: 

C2 = 2 log (^) [J(tmax) - J(Wn)] 



m" 



-1 



(13) 



(14) 



which goes to 1 in the typical case where y is small. The final integral is best 
performed in terms of the variable z = sy/w?: 

Sa\ f s \ log(l + z) 

and, taking the upper limit on z to infinity for simplicity, we finally have 

'1 + zq\ log(l + 2;o) 



a^ 



o-app = -^ log 



m? Km? 



log 



Zq / Zq 



(16) 



where zq is the lower limit on z. In practice, we want to impose a lower limit 
on fc° rather than on z, or y, since g^ '^iH usually be extremely close to its 
maximum. However, for fixed k^, the value of y is bounded by 

which is reached when p2 and k are coUinear and opposite to q2. So, unless 
K, is very close to 1, we have a lower bound 

z>zo = -^^^ , (18) 

and we shall use this in the event generator. Note, however, that upon gen- 
erating events with this lower bound, we unavoidably also get events with 
fc° < kE. These we discard by putting their weight to zero: even without 
additional experimental cuts, therefore, the generated event sample will con- 
tain a modest fraction of events with zero weight. 

We are now in a position to generate the events. This is done by ap- 
plying mappings in the order opposite to that of the integration, namely, 
first z, distributed according to Eq.(^), then t, according to Eg. ([T^) , and 
finally the angular variables, of which only v.y is not completely trivial. One 
minor remark is in order here: the generation algorithm assumes an infinite 
upper limit on z. Very occasionally, therefore, a y value will occur that is 
unphysically large. This is immediately indicated by tyj^[-^ becoming com- 
plex, in which case the event weight is put to zero and the event skipped. 
After having obtained all the phase space variables, we can construct the 
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momenta of the outgoing particles. This is in principle straightforward, but 
is made tricky by the enormous Lorentz boost necessary to go from the lab 
frame to the i?-frame and back. An ample collection of numerical pitfalls 
occurs, which we avoid by using expressions in which the major cancellations 
(from order 1 down to order m? / E'^) have been performed by hand. Finally, 
we calculate M and Mapp, also with an eye to the numerical stability. For 
the computation of M in an numerically stable manner, we refer to 0. In 
addition, since |t| is usually very small, t^ will occasionally lead to underflow 
problems. We therefore evaluate t^M and t^Mapp instead. Finally, we re- 
mark that ordinary double precision (REAL*8) suffices for our numerical work; 
extended precision (REAL* 16), such as employed in [§], is not necessary. 

3 Running BBBREM 

The Monte Carlo has a very limited number of input parameters, read at 
initialization from the standard input unit. They are, in order: 

1. ROOTS the total CMS energy of the incoming e^e~ system, in GeV. Any 
realistic value is allowed, from around 1 GeV (such as for the DA$NE 
machine) up to hundreds of GeV for NLC/CLIC. 

2. RKO the photon energy cutoff fraction k. In realistic cases, k will be 
of the order of percents, but it can be put much lower (although not 
to zero). When RKO increases, the fraction of internally rejected events 
(see the discussion above) increases slightly. RK0=1 is not allowed. 

3. NEVENT the requested number of generated events. This includes events 
that come out with zero weight. 

4. NRAN a flag determining the source of (pseudo-)random numbers used. 
In principle, any reliable such source may be substituted instead of the 
routine RANDOM, which occurs at 8 places in the generator code. We 
provide two alternatives: 

• NRAN=1 A very crude, low-period multiplicative congruential algo- 
rithm. Its only merit is portability, and it should only be used for 
checking the program against the test run mentioned below. 



• NRAN=2 A simple additive quasi-random number algorithm. We 
have used this to obtain the numerical results mentioned below. 
It is, however, not strictly portable since different rounding pro- 
cedures on different machines lead to different sequences. The 
results below should, however, be reproduced to within statistical 
accuracy. 

The input parameters ROOTS and RKO are transferred to the event generator 
subroutine BCUBE by common 

COMMON/EXPERC/ ROOTS, RKO 

and NRAN is transferred to the random number source RANDOM by common 

COMMON/RANCHN/ NRAN 

When called the first time, BCUBE initializes, and prints the result for (Tapp, 
in millibarn. Note that this is not intended as an accurate numerical approx- 
imation to the actual cross section! Upon each call, BCUBE fills the common 

COMMON/LABMOM/ Pl(4) ,P2(4) ,Q1(4) ,Q2(4) ,QK(4) ,D1 ,D2,T, WEIGHT 

which contains the essential information on the event: the four-momenta Pi, 
P2, Qi, q2 and k^^ (the fourth component being the energy), all in GeV; the 
invariants Ai, A2, and t, all in GeV^, whose computation, directly from the 
returned momenta, is numerically unstable; and, finally, the event weight, in 
millibarn. As given, the main program yields, after the generation of NEVENT 
events, the following significant numbers: 

1. RANCHK the result of one last call to RANDOM. This serves as a check 
on the portability of the program, since it depends both on the actual 
performance of the random number generator, and the number of ran- 
dom numbers that was actually used — this number is variable due to 
various decision and rejection steps. 

2. SMC the computed cross section, i.e. the average weight. 

3. SER the estimated error on the value of SMC, defined in the standard 
way, using the variance of the weight distribution. 
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In addition, some information on the weight distribution is given: Wl, the sum 
of the weights; W2, the sum of their squares; WNIL, the number of events with 
zero weight; WNEG, the number of events with negative weight: these should 
not occur, but if they do this would indicate rounding problems, or the use of 
an inappropriate form for M (for instance, negative weights arise if the form 
of ref. is used, although the cross section is still essentially correct); and, 
finally, WMAX, the maximum weight occurring in the generation. Concerning 
this last number it should be remarked that usually WMAX is quite a bit larger 
than SMC: the weight distribution has a low tail. This can be traced back 
to small values of Ci. In principle, the t integral in Eq.(|T2|) can also be 
performed exactly, without the benefit of Ci, leading to a slight complication 
in the algorithm. This might be desirable in such cases as necessitate the 
used of unweighted events, since high values of WMAX/SMC reduce the efficiency 
for unweighted events; for instance, if the produced momenta are to be used 
for tracking in machine studies. There, however, the tracking itself takes so 
much time that a somewhat inefficient unweighted-event generation is not a 
problem. We have therefore decided for the simpler algorithm. In any case, 
the high-weight tail does not prohibit the attainment of a small Monte Carlo 
error estimate. 



4 Results 

We shall now discuss some results obtained with BBBREM. In the first place, 
we reproduce here the output of a test run with the following choice of pa- 
rameters: R00TS=91.2d0, RKO=0.01dO, NEVENT=10000, and NRAN=1. It is: 



*** bbbrem : beam-beam bremsstrahlung *** 
*** authors r. kleiss and h. burkhardt *** 

total energy 91.200000 gev 

minimum photon energy 0.010000 * beam energy 

approximate cross section 0.627848D+03 millibarn 
random: very crude random numbers chosen 
******** cross section evaluation ******** 
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random number check 0.1471576690674 
total # of events 10000. 

sum of weights 0.323102D+07 

sum of (weight**2)s 0.390034D+10 

max weight occurred 0.221525D+05 

# of zero weights 1805. 

# of negative weights 0. 
computed cross section 0.323102D+03 millibarn 

+/- 0.534452D+01 millibarn 

We now turn to a more physical application of BBBREM, concerning the life- 
times of e+e~ colliding beams. As discussed in [|I|, the process ([I|) is actually 
the dominant process by which electrons and positrons are lost from the 
beams in storage rings such as LEP: an energy loss of an electron/positron 
due to bremsstrahlung that exceeds the so-called r.f. bucket half-height Sb 
cannot be made up for by the r.f. system, and the particle will be lost. At 
LEP and most other e~^e~ machines, Sb is typically about 10~^. One can 
therefore calculate, from first principles, the corresponding cross section: at 
^/s = 91.2 GeV, corresponding to LEP running at the Z^ peak, one finds 
then a ~ 320 mb, both from BBBREM and analytical calculations such as [^]. 
Actual measurement [jT| finds, however, a cross section of about 210 mb. As 
discussed in detail in [|l|, we ascribe this reduction to a density effect: in- 
side each bunch, the electromagnetic fields from all electrons (say) overlap in 
such a way that each electron has a finite interaction range, in contrast to the 
usual infinite range. This range is of the order of half the average distance 
d between two electrons in the bunch (at rest!), which is |]I| about 3.3 /im, 
corresponding to a momentum-transfer squared {hc/d)"^ = 3.56 ■ 10^^^ Gev^. 
Denoting this number by tc, we may then impose a cutoff on the momentum 
transfers allowed in the process (|ID, in two ways: 

• hard cutoff: events with |t| < tc are forbidden. In the Monte Carlo 
program this just means putting their weight to zero. 

• soft cutoff: the photon mediating the positron scattering is assigned 
an effective mass hc/d, leading to an effective exponentially suppressed 
electromagnetic potential of the electrons. In the Monte Carlo, this is 
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very simply iinplemented by multiplying the event weight with a factor 

Note that it would be very hard to implement either of these two cut-off 
procedures in any way but that of Monte Carlo. Below we present results for 
four different values of ROOTS: 1.02 GeV, corresponding to DA$NE at the 
$-meson peak, 60 GeV, appropriate to TRISTAN, 91.2 GeV, corresponding 
to LEP 1 at the Z" peak, and 190 GeV, typical for LEP200. In each case 
we have used RK0=0.01, NEVENT=1000000, and NRAN=2. The cross sections 
are given in millibarn, for each of the three cases: no cutoff, hard cutoff, soft 
cutoff. 



ROOTS (GeV) 


no cutoff 


hard cutoff 


soft cutoff 


1.02 


208.1 


198.7 


194.6 


60 


308.6 


210.2 


204.0 


91.2 


318.9 


210.1 


203.9 


190 


336.8 


209.8 


203.7 



The estimated error is 0.4 millibarn for DA$NE, and 0.5 millibarn in the 
other cases. We see that the LEP 1 cross sections with the cutoff are in 
much better agreement with the observation. For more details about this 
physical, rather than computational, issue we refer to M. 
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